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In the absence of synaptic coupling, two or more neural oscillators may become synchronized by 
virtue of the statistical correlations in their noisy input streams. Recent work has shown that the 
degree of correlation transfer from input currents to output spikes depends not only on intrinsic 
oscillator dynamics, but also depends on the length of the observation window over which the 
correlation is calculated. In this paper we use stochastic phase reduction and regular perturbations 
to derive the correlation of the total phase elapsed over long time scales, a quantity which provides 
a convenient proxy for the spike count correlation. Over short time scales, we derive the spike count 
correlation directly using straightforward probabilistic reasoning applied to the density of the phase 
difference. Our approximations show that output correlation scales with the autocorrelation of the 
phase resetting curve over long time scales. We also find a concise expression for the influence of 
the shape of the phase resetting curve on the initial slope of the output correlation over short time 
scales. These analytic results together with numerical simulations provide new intuitions for the 
recent counterintuitive finding that type I oscillators transfer correlations more faithfully than do 
type II over long time scales, while the reverse holds true for the better understood case of short 
time scales. 



While the jury is still out on the functional role of 
synchrony and correlations in neural firing, the ubiquity 
of these phenomena in the nervous system is suggestive. 
One long-standing hypothesis holds that correlated activ- 
ity in the visual system underlies feature binding. Syn- 
chronous oscillations may also play a role in amplifying 
signals [T], transmitting information from one layer to 
another [2Hl]i or such oscillations may encode informa- 
tion directly [51-fT^ . On the other hand, correlations may 
negatively impact the signal-to- noise ratio [TS'-T^, and 
excessive synchrony is a hallmark of neurological disor- 
ders such as epilepsy and Parkinson's disease. 

To understand the function of oscillatory correlations, 
or one day achieve clinically relevant control over them, 
we must first understand the underlying biophysical 
mechanisms. While synchrony can arise as the result of 
anatomical connectivity between neurons, much recent 
work [171 - 122] has brought to light ways in which corre- 
lated activity develops from the inherent stochastisicity 
of neural systems. Thus, in the absence of direct cou- 
pling, two or more neural oscillators may become syn- 
chronized by virtue of the statistical correlations in their 
noisy input streams - a phenomenon we will refer to as 
stochastic synchrony. 

For our analysis of stochastic synchrony, we appeal to 
the theory of weak coupling, which holds in the stochastic 
context provided the amplitude of the noise is sufficiently 
small. In particular, a number of groups [TSHiniUS] have 
proved that the phase reduction technique |24) can be 
applied to oscillators receiving additive noise. Thus, we 
reduce a noisily driven oscillator to a scalar differential 
equation describing the evolution of the phase. This so- 
called phase equation depends only on the properties of 
the noise and the oscillator's phase resetting curve (PRC) 
which characterizes how small perturbations influence 
the oscillator's subsequent timing or phase. 

Neural oscillators can be classified into two types ac- 



cording to the bifurcations that occur as the dynamical 
system goes from a stable rest state to a stable limit 
cycle. Furthermore, the oscillator's bifurcation class has 
been shown to determine the shape of it's PRC and there- 
fore it's ability to synchronize. Type I oscillators undergo 
the saddlc-nodc-on-an-invariant-circle, or SNIC, bifurca- 
tion and the resulting PRC is strictly positive, indicat- 
ing that perturbations can only advance the oscillator's 
phase. Type II cells undergo the Andronov-Hopf bifur- 
cation, which produces a PRC with both negative and 
positive regions; typically, inputs occurring early in the 
cycle can delay the phase while later inputs advance it. 
See Fig.^. 

An expanding body of work has demonstrated that 
over short time scales of less than one period, type II 
oscillators are more susceptible to stochastic synchrony 
than type I. This has been shown via simulations and 
in vivo [171 125| , by deriving the probability distribution 
of the phase difference |26| , by minimizing the Lyapunov 
exponent of the phase difference [37], and most recently 
by calculating the spike count correlation over a range 
of time windows 28]. The latter study further reports 
that this finding reverses over long timescales, namely 
that type I oscillators transmit correlations more faith- 
fully than type II when observed over lengths of time 
much greater than one period. 

In Section |l] we provide a brief introduction to the 
phase reduction technique in a stochastic setting. Next 
in Section [IT] we use regular perturbations to give a novel 
and straightforward analysis of correlation transfer over 
long time scales. To facilitate our derivation, we use the 
total elapsed phase as a proxy for the spike count. Note 
that the total phase (modulo the period) and the spike 
count differ by at most one, which is a negligible quantity 
when many spikes have been observed over a long time 
window. The expression we derive for the correlation co- 
efficient of the total phase agrees both qualitatively and 
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FIG. 1. We use the parametrization A(^) = — sin(6' + a) + 
sin(a) to vary the PRC smoothly from type I (red), where a = 
^ and A{6) = 1 — cos{9), to type II (blue), where q = and 
A{0) — — sin(S). Note that intermediate values of a produce 
PRC shapes (dashed purple) that more closely resemble those 
found empirically in vivo. 



quantitatively with the results found in [2F. 

In Section IIIII we consider short time scales less than 
or equal to the period of the oscillation. In this case, 
the total phase cannot be used to approximate the spike 
count. We therefore derive the spike count correlation 
directly, using simple probabilistic reasoning applied to 
the density of the phase difference. Our analytic results 
together with Monte Carlo simulations corroborate ear- 
lier work showing type II oscillators transfer correlations 
more readily than type I over short time windows. 

I. NOISY OSCILLATORS 

Let us begin with a neural oscillator receiving additive 
noise with equations of motion given by 

dX = F{X)dt + a^, 

where X e M" and ^ is a white noise process. When a = 
0, we assume the noiseless system has an asymptotically 
stable periodic solution Xo{t) — Xo{t + T) with period t. 

As in the deterministic case, we can reduce this high- 
dimensional system to a scalar equation for the evolution 
of the phase 9 around the limit cycle. Let </> : K" — > §^ 
map a neighborhood of the limit cycle to the phase on 
a circle. That is, 6 = 0(X), with 9 e [0,1). Then 9 
satisfies 

^ = l + aVx</'(X).e, 

where we have normalized the unperturbed period to be 
one. Next we can close the equation by assuming the 
noise amplitude a is sufficiently small, so that the sys- 
tem trajectory can be approximated by the noiseless limit 
cycle Xq: 

9^1 + aZ{9)-^, (1) 

where Z{9) = V x4'{^o{(^)) is the adjoint, or phase- 
dependent sensitivity of the trajectory to perturbation 



along the limit cycle. In the case of a neural oscillator, 
we assume the noisy perturbations arise as the result of 
stochastic synaptic input, which influences only the volt- 
age variable. Hence Z(9) has only one nonzero compo- 
nent, which is proportional to the phase resetting curve 
Ai9). 

Thus far, wc have used the conventional change of vari- 
ables to obtain Eq.([l]), which therefore must be under- 
stood as a stochastic differential equation (SDE) in the 
Stratonovich sense. In order to eliminate the correlation 
between 9 and ^ we must use the Ito change of variables, 
which will introduce an additional drift term: 

2 

9 = l + aAi9)^+^A'i9)A{9). 

Here ' denotes differentiation with respect to 9. For a 
detailed discussion of phase reduction in noisy oscillators 
see ESI. 



II. CORRELATION TRANSFER OVER LONG 
TIME SCALES 

We now consider the transfer of correlations over time 
scales much larger than the natural period of the oscil- 
lators. Given the level of correlation between the noisy 
inputs, we wish to know what level of correlation remains 
between the spike count of two oscillators after some 
time. For analytic convenience, however, we will use 
the total phase that has elapsed as a proxy for the spike 
count. Since these quantities differ by at most one, the 
discrepancy will be negligible for the large spike counts 
that accrue over long time scales. 

Our system will consist of two identical phase oscilla- 
tors receiving weak, correlated, but not identical, addi- 
tive white noise. Keeping only terms up to order a, we 
have 

0'i = i + (7A(0i)a(<) 

02 = 1 + 'tA(02)6W- (2) 
The noise takes the form 

6 - Vcec + \/r^es, (3) 

where ^a, and are mutually independent, zero 
mean white noise processes, and c G [0, 1] is the cor- 
relation between and ^2, which we will refer to as the 
input correlation. 

Next let us rewrite Eq.([2]) in the form of integral equa- 
tions: 

6»i (t) = i + 9i {Q) + a f A{9i (s))^ {s)ds 

92(t) = t + 92{0) + a f A(02(s))6(s)ds. 
Jo 
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FIG. 2. The steady state distribution P{<j}) of phase differences </!> is shown for type I (red) and type II (blue) as well as for 
intermediate PRCs (dashed purple). Note that the unperturbed period of the oscillators is 2it. (A) Input correlation c = 0.4. 
(B) Input correlation c — 0.8. 



Let T be length of the window of time over which 
we will observe the system. Throughout this discus- 
sion we will assume that our system has reached equi- 
librium, and that time has been reparametrized so that 
our observation takes place on the interval t E [0,T]. 
In order to quantify the total phase traversed during 
this time, we subtract the initial phases by defining 
qi{T) = e^iT) - 6,{0) for i = 1,2. Thus the total phase 
traversed over a time window of length T is given by: 



q,{T)=T + a A((?,(s))e,(s)ds. 



with qi{0) — for i — 1,2. Finally, since we assume a 
is small, let us simplify the integrands by expanding the 
phase to lowest order: 



Then we have A (6*^(5)) = A{t + 6',;(0)), and thus 



q,{T) =T + a A{s + 0,(O))^,(s)ds 



(4) 



(5) 



When taking expectations of the quantities in Eq.([5]), 
we must keep in mind that there are four stochastic vari- 
ables over which averaging must take place. In partic- 
ular, we must average over the white noise signals ^i(t) 
and <^2(^) and the initial conditions Oi{Q) and 02(0) . 

Assuming we begin observation after the system has 
reached equilibrium, we can take one of the initial con- 



ditions, say 9i{0), to be distributed uniformly on the in- 
terval [0, 27r]. However, at equilibrium the phases obey 
the steady state probability distribution P{(t>) derived in 
p6] and [30] , which depends only on the phase difference 
(j){t) = 02{t) - 0i{t). Therefore, the average of Eq.^ is 
computed as 



1 

2^ 



T + a I A{s + x)^i{s)ds 




P{y-x)x 



= T 



T + a I A{s + x) ids)) ds 
P{y -x)x 





2tt /.2ir 



dxdy 







a 

T 

A{e,{s)) {C,{s)} dsdxdy 











(6) 



where 2tt is the unperturbed period of the oscillators, 
P(0) is the steady state probability distribution of the 
phase difference, and x and y represent the initial condi- 
tions 9i{0) and 02{O), respectively. The last line follows 
because the white noises have zero mean. 

Our goal is to compute the correlation of the total 
phase traversed by the two oscillators: 



Cor[gfi,g2] 



Cov[gi,g2] 
^Var[gi]Var[g2] 



(7) 
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First, we derive the covariance as follows 



Cov[qi,q2]{T) = E[(gi(r) - E[<7i(T)])(g2(r) - EfelT))]] 
= E[(gi(r)-T)fe(T)-r)] 

i-T rT 



E 



a- / A{s + ei{0))^i{s)ds A{s' + e2{o))Us')ds' 
Jo 

T rT 



'2 ™ / / P{y-x) / / A(s + a;)A(s' + 2/)(5(s-s')dsds'da;d?/ 
Iq Jo Jo Jo 

2 '-'in 



27r 



27r 



./o 



P{y — x) / A(s + cc) A(s + y)(is(ia;(i2/. 







Similarly, we find the variance to be 



Var[qi](r) = E[(gi(T)-E[gi(T)])2] 

P{y — x) I A(s + x)^dsdxdy. 



1 



27T 



Note that we therefore have Var[gi] — Var[(72], and 
hence the denominator of Eq.([7]) can be simplified: 
■yVar[(7i]Var[g2] — Var[gi]. This gives the total phase 
correlation as 

Cor[gi,g2](r) 

_ Jo'' Jo'' P{y - x) A(s + x)A{s + y)dsdxdy 



JT Jo'' - •^) Jo + xYdsdxdy 



(8) 



Now let h{x) — Jo'' ^{y)^{y+x)dy be the autocorrela- 
tion of the PRC, and let ~ 02{t) — 9i{t) represent the 
phase difference as before. Then we can rewrite Eq.(l8|) 



Cor[gi,g2](T) -c 



j^^p{mo)d^ 



Note that the right hand side no longer depends on T af- 
ter we switched the order of integration and canceled the 
resulting factors of T in both numerator and denomina- 
tor. Next we can do away with the denominator entirely, 
since h{0) does not depend on cf), which leaves simply 



(9) 



An expression for the steady-state probability density 
of the phase difference P{x) was derived by Marella and 
Ermentrout in (26j. Specifically, we have 
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FIG. 3. Output correlation for large time windows is shown 
as a function of the PRC shape parameter a. Note that when 
a = the PRC is a pure sinusoid and therefore the osci llat or 
is type II; when a = it/2, the oscillator is type I (see Eq.( 10 1). 



Theoretical curves (solid) are a good match for both the simu- 
lated total phase correlation (dotted) and the simulated spike 
count correlation (starred). Colors indicate the level of input 
correlation: 0.2 (blue), 0.4 (green), 0.6 (red), 0.8 (cyan), 0.99 
(purple). In all cases, noise amplitude a = 0.05. 



where G{x) = 1 — c {h{x)/h{0)), and iV is a normalizing 
constant, iV = 1/ J^^ 1/G{x)dx. Let us further define 
the PRC to be 



A(^; a) — — sin{9 + a) — sin(a). 



(10) 



where a is a parameter that allows us to vary the PRC 
shape smoothly between type I (a — 7r/2) and type II 
{a — 0). See Fig.([T]). Using this, the phase distribu- 
tion over long time scales becomes a function of input 
correlation and the PRC shape parameter: 



5 




a 



FIG. 4. The perturbation expansion of Cout for small input 
correlation (dashed) agrees well with the full output corre- 
lation (solid). Note that, to lowest order in dn, the output 
correlation goes to zero as the PRC shape parameter a goes 
to zero, that is, as the PRC shape approaches the pure type 
II. Colors indicate the level of input correlation: 0.01 (blue), 
0.05 (green), 0.1 (red). 



P{(j>; c,a) 



y/{c - l)(cos(2a) - 2)(2 + (c - 1) cos(2a)) 



27r(2 - c + (c - 1) cos(2a) - ccos((/»)) 



(11) 



In the special cases where a = 7r/2 and a = 0, Eq.( 10 ) 
and Eq.(ll), together with Eq.(|8|, yield 



Type I 

Ai{x) = 1 — cos(a;) 

, 73 Vc2 -4c + 3 
Pi {(f); c) = 



27r (3 - 2c- ccos((/>)) 



(12) 



Cout,I = 1 - ^ a/3(c- 3)(c- 1) 



Type II 

Aii{x) = — sin(x) 
1 



Pii{<i>;c) = 



2n (1 — ccos((/))) 



(13) 
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Our intuition for this finding can be honed by per- 
forming a further perturbation expansion, now assuming 
small input correlation. For sufficiently small c, we can 
make the approximation 



Gix) 



1 - c 



h(x) 
h{0} 



1+C 



h{x) 



When we substitute this into Eq.([9]) we find 



N 



^out 



2tt 



(14) 



where iV = 1/ ^"^^ {1 + Cinh{x) / h{Q)) dx is likewise ap- 
proximated to lowest order in c. 

The form of Eq. ([m]) demonstrates that output correla- 
tion scales with the integral of the PRC autocorrelation. 



and for the parametrized PRC in Eq.(10l we have 



2-K 



h{(j))dcj) = 47r2 sin(a)^ 



In particular, a = for the type II PRC, and hence 
Cout = to lowest order. Clearly, we have nonzero au- 
tocorrelation for nonzero a < f , and hence PRCs that 
deviate from pure type II will produce higher output cor- 
relation over the long timescales considered here. 

Expanding the remaining terms in Eq.( 14 ), we find the 



approximated output correlation takes the form 



Cout 



2csin(a)^ 



2 + c - (1 c) cos(2a) ' 



(15) 



In Fig.(|4]) we see that this approximation agrees with 
Eq.(|8]) for c = 0.01 and 0.05 but diverges for c = 0.1. 
Note that these curves would all lie below the lowest 
curve plotted in Fig.([3| if shown on the same scale. 

We verify the preceding analysis by simulating two 
phase oscillators perturbed by additive white noise as 
described in Eq.([2]) and Eq.(3|. The simulations used 
noise amplitude a = 0.05, and the input correlation took 
the values c e {0.2, 0.4, 0.6, 0.8, 0.99}. 

We computed the correlation coefficient of both the 
total phase and the spike count, using a range of obser- 
vation windows T. As shown in Fig.(3|, the total phase 
correlation and the spike count correlation agree closely 
both with each other and with the theoretical curves as 
a function of the PRC shape parameter a. 



As in [55], we see in Fig.([3]) that type I oscillators dis- 
play greater output correlation than type II oscillators 
for any fixed value of the input correlation c, a surprising 
finding in light of earlier results that demonstrated the 
opposite relationship over short windows of observation. 



III. SHORT TIME SCALES 

Now we will calculate the spike count correlation di- 
rectly for observation windows T that are shorter than 
or equal to the natural period, which we will assume to 
be 27r. First let us consider the probability that a spike 
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FIG. 5. Joint spiking probability for two oscillators receiving 
partially correlated noise is shown for observations windows 
T < 1, where 1 is the natural frequency of the oscillation. 
The subscripts ij indicate the probability that oscillator i or 
j does (1) or does not (0) spike. 



occurs in [0,T]. We say that oscillator i spikes when its 
phase 9i reaches 27r, which is to say Oi{T) > 2tt. As- 
suming as usual that the noise amplitude a is small, we 
expand the phase to lowest order as in Eq.Q, that is 
6»,(T) = e^{0)+T + O{a). Therefore the probability that 
oscillator i spikes is simply 



P[9, spikes] ^P[e,+T> 27r] 
P[e, does not spike] = P[6', + T < 27r]. 



For two oscillators, there are four possibilities for the 
joint spike count: 



P[9i does not spike, 6*2 does not spike] 

= P[6'i+r<27r,6l2+T<27r] 
P[9i spikes, 6*2 does not spike] 

= P[9i +T >2tt,92+T < 2n] 
P[9i does not spike, 6*2 spikes] 

= P[6'i + r < 27r,6»2 + T > 2tt] 
P[9i spikes, 6'2 spikes] 

= P[e'i +T >2tt,92+T > 2n]. 



These probabilities can be obtained directly by in- 
tegrating the density of the phase difference, Eq.(ll), 
over the appropriate domain. Note that this gives four 
discrete joint probabilities for each observation window 
T € [0,27r]. For convenience, let us define the following 



functions of T: 
fooiT) : 

foiiT) : 

fwiT) : 



P[Oi < 27r-r,e'2 < 2tt-T] 

2TT-T r2TT-T 

P{y — x)dxdy 



1 

27r JO JO 
P[6li > 27r-r,6'2 < 2tt-T] 

It: 

/ P{y — x)dxdy 

2iT-T JO 

P[9i <2tt-T,92> 2tt~T] 

27r-T /.2ir 

P{y — x)dxdy 



1 

2^ 



1 

Jo J2TT-T 
hi{T) := P[9i > 27r - T, 02 > 27r - T 

1 



2ti 



2-K-T J27r-T 



P{ij — x)dxdy. 



Let X be the random variable such that X = 1 if 6*1 
spikes during the observation period T, and X ~ ii 9i 
does not spike. Similarly, let Y represent the presence or 
absence of a spike in oscillator 92- Then the covariance 
is given by Cov[X, Y] = E[XY] - E[X]E[Y]. In terms of 
the functions defined above we have 

E[X] =0-(/oo + /oi) + l-(/io + /ii) 
= (,/io + /ii)=E[X2] 

Em =0-(/oo + /io) + l-(/oi+/ii) 
= (/oi + /ii)=E[r2] 
E[XY] - • • /oo + 1 • • /lo + • 1 • /oi + 1 • 1 • /n 
= /ii- 

A few simplifications are possible. In particular, the 
sum /io(T) + fii(T) is just the marginal probability that 
9i spikes within time T. Since 9i is uniformly distributed, 
this probability is simply Furthermore, we also have 
fi o = /oi by the symmetry of the density P, and hence 
■y/Var[A']Var[F] = Var[X]. Therefore the spike count 
correlation over short time windows is 



Coy[X, Y]{T;c) 
_ E[XY] - E[X]E[y] 
~ Var[X] 

_ fll — ifw + /ll)^ 

" (/io + /ii)(l-(/io + /ii)) 



(16) 



27r \^ 27tJ 



1 



2t:T - T2 



2n / P{y~x)dxdy-T^ 

J27r-TJ27r-T 



(17) 



Fig.([6jA_,B) shows how this analytically derived output 
correlation compares with numerical simulations for type 
I and type II oscillators, respectively. 
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FIG. 6. (A,B) Theoretical (solid) and simulated (dotted) output correlation curves are shown as a function of the observation 
window T < 27r. (A) Type I oscillators. (B) Type II oscillators. (C,D) The initial slope (dashed) of the spike count correlation 
(solid) is the linear approximation of Eq.([T7| at T = 0, which is given in Eq.( |19[ ). (C) Type I oscillators. (D) Type II oscillators. 
For all plots, noise amplitude a = 0.04, and colors indicate the level of input correlation; 0.2 (blue), 0.4 (green), 0.6 (red), 0.8 
(cyan), 0.99 (purple). In all cases, noise amplitude a = 0.05. 



We can make a further simplification by considering 
the linear part of Eq.(17l for T close to zero: 



Cout = T P(0) 



2ti 



Cout J I 




(18) 
(19) 



Thus, the initial slope of the output correlation is pro- 
portional to the peak of the stationary phase difference 
distribution, P((/))|0=o- Substituting P/(0) and Pii{0) 



from Eq.(12| and Eq.(13l, we obtain 



From here, it is clear that the initial slope of Cout is 
greater for type II than for type I oscillators; in fact the 
type II output correlation rises three times faster than 
the type I, to lowest order in c. See Fig.([6]G,D). 
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FIG. 7. Output correlation is shown as a function of 
intermediate-length observation windows T. Colors indicate 
the level of input correlation: 0.2 (blue), 0.4 (green), 0.6 (red), 
0.8 (cyan), 0.99 (purple). (A) Type II oscillators (solid) ex- 
hibit higher output correlations over short time scales than 
do type I (dashed). (B) This result reverses over long time 
scales. In all cases, noise amplitude a — 0.2. 



DISCUSSION 

We have demonstrated a novel approach to approxi- 
mating the spike count correlation of noisy neural oscil- 
lators over both long and short time scales. In the case 
of long windows of observation T much greater than the 
natural period of oscillation, we used the total elapsed 
phase, modulo the period, as a proxy for the spike count. 
The difference between these quantities is at most one 
and hence is negligible for when many spikes are ob- 
served over large time windows T. In our perturbation 
expansion to lowest order in the noise amplitude, a, the 
correlation between oscillators depends only on the PRC 
and the stationary distribution of the phase difference. A 
further approximation assuming small input correlation 
c reveals that output correlation scales with the autocor- 
relation of the PRC, which is a nonnegative quantity that 
equals zero precisely when the PRC is a pure sinusoid, 
i.e., when the oscillator displays type II dynamics. This 
observation sheds some light on the counterintuitive find- 
ing, first reported by Barreiro, et al. whereby type 
I oscillators transfer correlations more faithfully than do 
type II over long time scales, although the reverse holds 
true for the better understood case of short time scales. 

Using straightforward probabilistic reasoning, we com- 
puted the spike count correlation directly for short time 
scales. In the limit of small T and small c, we obtain 
an expression for the initial slope of the output correla- 
tion, also known as the correlation susceptibility [8 . In 
[8], de la Rocha, et al. use a phenomenological model to 
explore the complex relationship between susceptibility, 
firing rate and threshold nonlinearities. The present anal- 
ysis illustrates the contribution of bifurcation structure 
via phase resetting dynamics. In particular, the suscepti- 
bility is proportional to the peak of the stationary phase 
difference distribution, P(0)|0=O; which in turn depends 
on the shape of the PRC. 

Our analytic expressions in the limit of small noise 
agree well with spike count correlations computed from 
simulated oscillators. However, for tractability we in- 
cluded only terms of order one in the perturbation ex- 
pansion of the phase given in Eq.(|4]). As a result, the 
present analysis cannot account for the slow drift of the 
correlation due to noise, which is visible for values of T 
near 2n in Fig. ([6]), and is even more apparent for the in- 
termediate values of T shown in Fig.([7]). New analytic 
methods capable of addressing non-extremal cases would 
shed light on this and many other questions in mathe- 
matical biology. 
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